Leakage Belt Delineation

VM0048-compliant delineation of valid leakage area for the Gola REDD+ Forest Carbon Project project area expansion

Seamus Murphy mailto:Seamus.Murphy@winrock.org (Winrock International)https://winrock.org/ , Krysla Lima Fernandes mailto:Krysla.Lima@winrock.org (Winrock International)https://winrock.org/ , Rajesh Bista Rajesh.Bista@winrock.org (Winrock International)https://winrock.org/ , Meyru Bhanti Meyru.Bhanti@winrock.org (Winrock International)https://winrock.org/
2025-05-17

Project Area & Jurisdictions

country = sf::read_sf("./data/AOI/liberia_boundary_national.shp") |>sf::st_transform(32629)
counties = sf::st_read("./data/AOI/places_poly_county.shp") |>sf::st_transform(32629)
Reading layer `places_poly_county' from data source 
  `/Users/seamus/repos/rspb-redd-risk-new/data/AOI/places_poly_county.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 16 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -11.50675 ymin: 4.353908 xmax: -7.367323 ymax: 8.551925
Geodetic CRS:  WGS 84
jurisdiction = counties |>dplyr::filter(name=="Grand Cape Mount County"|name=="Gharpolu County")
jurisdiction$name = 'Grand Cape Mount & Gharpolu Counties'

aoi = sf::read_sf("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/AOI/Archive/ProjectArea_CF-Expansion_051525/updated_areas.shp") |>
  sf::st_make_valid() |>
  sf::st_transform("EPSG:32629") |>  
  sf::st_cast("MULTIPOLYGON") |> sf::st_as_sf()# |> dplyr::select("Name")

aoi$area_ha = round(as.numeric(sf::st_area(aoi) * 0.0001, 4))
aoi |> sf::st_drop_geometry() |> janitor::adorn_totals() |> 
  flextable::flextable() |> 
  flextable::fontsize(size=8,part="all") |> 
  flextable::autofit()

Name

Shape_Leng

Shape_Area

area_ha

Lukasu

138,549.68

469,973,056

46,997

Mbarma CF

91,209.87

443,230,136

44,323

Lower Sokpo CF

59,692.55

140,691,204

14,069

Upper Sokpo CF

62,831.33

109,225,715

10,923

Zuie CF

92,221.99

362,360,113

36,236

Gola Forest National Park

182,951.78

890,423,315

89,042

Tonglay CF

76,824.36

195,953,085

19,595

Norman CF

77,870.76

123,489,968

12,349

Foya

351,458.82

1,048,825,566

104,885

Kposo

77,241.08

307,632,539

30,763

Gola Konneh CF

117,719.73

531,018,504

53,102

Gola Konneh CF (ext. S)

56,323.80

108,968,871

10,897

Gola Konneh CF (ext. N)

96,399.13

285,898,112

28,590

Lower Sokpo CF (ext.)

45,714.33

47,712,283

4,771

Maima CF (ext.)

59,344.81

139,163,854

13,916

Zuie CF (ext.)

17,920.40

14,953,753

1,495

Yangaya

46,462.77

58,336,055

5,834

Koninga

71,795.27

317,410,805

31,741

Lobarsu

62,061.63

165,254,549

16,548

Bade

99,147.36

513,666,252

51,367

No-mans land

209,804.70

394,466,117

39,447

Lower Guma

56,222.96

142,192,367

14,219

Central Guma

31,669.60

38,501,602

3,850

Upper Guma

54,605.75

137,757,843

13,771

Hassala

62,847.92

145,701,577

14,570

Buluyeama

47,779.35

100,212,025

10,021

Hembeh

199,395.60

461,835,431

46,184

Total

2,546,067.34

7,694,854,697

769,505

Derive Leakage Area Belt

aoi_union = sf::st_transform(aoi, 32629) |>  sf::st_union() |> sf::st_make_valid()
leakage_buffer = sf::st_buffer(aoi_union, dist = 5500, endCapStyle="ROUND") |>
  sf::st_as_sf()
leakage_buffer = concaveman::concaveman(leakage_buffer, concavity=5) |> 
  sf::st_zm() |> sf::st_make_valid()
leakage_belt_whole = sf::st_buffer(leakage_buffer, dist = 4500, endCapStyle="ROUND") |>
  sf::st_make_valid() 

tmap::tmap_mode("view")
tmap::tm_shape(leakage_belt_whole) + 
  tmap::tm_polygons(col="orange",fill="orange",fill_alpha=0.3,lwd=1)+
  #tmap::tm_add_legend(type="lines",col="orange",labels="Leakage Belt (10km)") +
  tmap::tm_shape(aoi) + tmap::tm_borders(lwd=1, col="white") +
  tmap::tm_basemap("Esri.WorldImagery")
Reading layer `LeakageBelt_10k-Radius_UnFiltered-Whole' from data source `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered-Whole.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -11.31102 ymin: 7.024169 xmax: -9.847199 ymax: 8.24916
Geodetic CRS:  WGS 84
leakage_belt    = sf::st_difference(leakage_belt_whole, st_union(st_combine(aoi_union)))
leakage_belt$area_ha = round(as.numeric(sf::st_area(leakage_belt) * 0.0001, 4)) 
leakage_belt = sf::st_intersection(country, leakage_belt)

tmap::tmap_mode("view")
tmap::tm_shape(leakage_belt) + 
  tmap::tm_polygons(col="yellow",fill="yellow",fill_alpha=0.3)+
  #tmap::tm_add_legend(type="lines",col="yellow",labels="Leakage Belt (10km)") +  
  tmap::tm_shape(aoi) + tmap::tm_borders(lwd=0.4, col="white") + 
  tmap::tm_text(text="Name", size=0.5, col="white") +
  tmap::tm_basemap("Esri.WorldImagery")

# save locally
#sf::wt_write(leakage_belt, "OneDrive.../20087 - RSPB Gola Feasibility/Deliverables/
#  Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered.zip") 
#sf::st_write(leakage_belt_whole, "OneDrive.../20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered-Whole.shp")
Reading layer `leakage_belt_10km_liberia_ext' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered/leakage_belt_10km_liberia_ext.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 245117.9 ymin: 776808.3 xmax: 406590.1 ymax: 908414.4
Projected CRS: WGS 84 / UTM zone 29N
Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 245117.9 ymin: 776808.3 xmax: 406590.1 ymax: 908414.4
Projected CRS: WGS 84 / UTM zone 29N
                      Name Shape_Leng Shape_Area
1 Leakage Belt 10km Radius   138549.7  469973056
                        geometry area_ha
1 MULTIPOLYGON (((245306.7 79...  298838

Derive Leakage Area Masks

Roads Mask

roads_one = sf::st_read("~/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_one.shp")
roads_two = sf::st_read("~/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_two.shp")

# we have simplify mask shapefiles and split them up to shorten computing 
# time & avoid crashing. See option for "harsh" simplification on line 163
roads_one_simplified = roads_one |> sf::st_make_valid() |> sf::st_cast("MULTILINESTRING") |> 
  rmapshaper::ms_simplify(keep=0.5)
roads_two_simplified = roads_two |> sf::st_make_valid() |> sf::st_cast("MULTILINESTRING") |> 
  rmapshaper::ms_simplify(keep=0.5)

# bigger file needs more simplificaiotn
roads_one_simplified_harsh = rmapshaper::ms_simplify(
  roads_one_simplified, keep=0.01) 

# now apply buffer operation, but note this takes time. Its 
# advised processing inputs as much as possible before running
roads_one_buffer = sf::st_buffer(
  roads_one_simplified_harsh, 
  dist = 10000, 
  nQuadSegs = 5,
  endCapStyle="ROUND", 
  joinStyle = "ROUND",
  mitreLimit = 1,
  singleSide = FALSE
  )

roads_two_buffer = sf::st_buffer(
  roads_two_simplified, 
  dist = 10000, 
  nQuadSegs = 5,
  endCapStyle="ROUND", 
  joinStyle = "ROUND",
  mitreLimit = 1,
  singleSide = FALSE
  )

# Combine, dissolve and cast to single feature
roads_mask = sf::st_combine(roads_one_buffer, roads_two_buffer) |>
  sf::st_union() |> sf::st_cast("POLYGON")

# Visual check
tmap::tmap_mode("view")
tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=0) +
  tmap::tm_shape(roads_one_simplified_harsh) + tmap::tm_lines(lwd=2, col="red") +
  tmap::tm_shape(roads_two_simplified) + tmap::tm_lines(lwd=2, col="green") +
  tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=1, col="pink") + 
  #tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_scale_bar(position = c("RIGHT", "BOTTOM"), text.size = .5) + 
  tmap::tm_compass(color.dark = "gray60", text.color = "gray60", position = c("left", "top"))

# Save output to MASKS folder and purge memory
#sf::st_write(roads_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASKS/LeakageMask_Roads_10km-Buffer_051625.shp", delete_dsn=T)
Reading layer `roads_simplified_one' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_one.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 842 features and 16 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 237648.4 ymin: 761962.5 xmax: 400449.6 ymax: 923370.5
Projected CRS: WGS 84 / UTM zone 29N
Reading layer `roads_simplified_two' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_two.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 255 features and 16 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 237648.4 ymin: 762265 xmax: 400334.4 ymax: 923370.5
Projected CRS: WGS 84 / UTM zone 29N
Reading layer `Road_Mask_10km-Proximity_051625' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/Archive/LeakageMask_Roads_10km-Buffer_051625/Road_Mask_10km-Proximity_051625.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 17 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 227660.4 ymin: 752279.5 xmax: 410286 ymax: 933257.1
Projected CRS: WGS 84 / UTM zone 29N

Habitat Mask

# import inputs
wetlands = terra::rast("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Wetlands/GLWD_EPSG32629.tif")
protected_areas = sf::st_read("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Protected Areas/Archive/WDPA_Mar2025_Public_32629_GOLA.shp")
Reading layer `WDPA_Mar2025_Public_32629_GOLA' from data source 
  `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Protected Areas/Archive/WDPA_Mar2025_Public_32629_GOLA.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 6 features and 30 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 240189.7 ymin: 751293.6 xmax: 411102.7 ymax: 895314.1
Projected CRS: WGS 84 / UTM zone 29N
leakage_belt_crop = sf::st_transform(sf::st_as_sf(leakage_belt_whole), 32629) |> terra::vect()
wetlands          = terra::crop(wetlands, leakage_belt_crop, mask=T)

# tidy labeling
code_dict_2 <- data.frame(
  id = c(1, 4, 7, 10, 12, 14, 15, 18, 20, 21, 26, 31),
  label = c(
    "Freshwater lake",                              # 1
    "Large river",                                  # 4
    "Small streams",                                # 7
    "Riverine, regularly flooded, forested",        # 10
    "Riverine, seasonally flooded, forested",       # 12
    "Riverine, seasonally saturated, forested",     # 14
    "Riverine, seasonally saturated, non-forested", # 15
    "Palustrine, seasonally saturated, forested",   # 18
    "Ephemeral, forested",                          # 20
    "Ephemeral, non-forested",                      # 21
    "Tropical peatland, forested",                  # 26
    "Other coastal wetland"                         # 31
  ))

levels(wetlands) <- code_dict_2
wetlands[wetlands == 0] <- NA

# derive wetland mask
wetlands_mask <- wetlands
wetland_classes <- c(1, 4, 7, 10, 12, 14, 15, 18, 20, 21, 26, 31)
terra::values(wetlands_mask) <- ifelse(terra::values(wetlands) %in% wetland_classes, 1, NA)

# save locally for faster computing
#raster::writeRaster(wetlands_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_Wetland-GLWD_051625.tif", overwrite=T)
#sf::st_write(protected_areas, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_ProtectedAreas_WDPA_051625.shp", delete_dsn=T)

Slope Mask

# skipping these operations here (est. time 12 mins)
DEM = terra::rast("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/DEM/DEM_EPSG32629.tif") 

# derive slope percentage from degree 
slope_degrees = terra::terrain(DEM, v="slope", unit="degrees")
slope_percent = tan(slope_degrees * (pi / 180)) * 100
slope_percent = terra::clamp(slope_percent, 0, 100) 
slope_invalid = slope_percent > 10
slope_invalid[slope_invalid == 0] <- NA
slope_mask = terra::as.polygons(slope_invalid, dissolve=T)|>sf::st_as_sf()|>sf::st_union()

# save locally & reload
sf::st_write(slope_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_Slope10%-Invalid_051625.zip", delete_dsn=T)
slope_mask = sf::st_read("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/Archive/LeakageMask_Slope10%-Invalid_051625/slope_poly_simplified.shp")

Visual Check

# Visual check
tmap::tmap_mode("view")
tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=0) +
  tmap::tm_shape(wetlands) + tmap::tm_raster(col.legend = tm_legend("Wetlands (GLWD")) +
  tmap::tm_shape(leakage_belt) + tmap::tm_polygons(col="yellow",fill="yellow",fill_alpha=0.4, lwd=1.5)+ 
  tmap::tm_shape(protected_areas) + tmap::tm_polygons(fill="ORIG_NAME", fill.legend = tm_legend("Protected Areas (WDPA)")) +
  tmap::tm_shape(aoi) + tmap::tm_borders(lwd=1.5, col="red") + 
  tmap::tm_text(text="Name", size=0.3, col="black") +
  tmap::tm_shape(roads_mask) + tmap::tm_borders(col="purple") +
  #tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = .5) + 
  tmap::tm_compass(color.dark = "gray60", text.color = "gray60", position = c("left", "top"))

Apply Leakage Area Masks

# clip
roads_leakage_mask = sf::st_intersection(leakage_belt, roads_mask)
slope_leakage_mask = sf::st_intersection(leakage_belt, slope_mask)
wetlands_leakage_mask = sf::st_intersection(leakage_belt, wetlands_mask)
protected_areas_leakage_mask = sf::st_intersection(leakage_belt, protected_areas)

# merge. Btw actually makes better sense to keep these seperate. 
# They are easier to operate seperately due to their size and linework.
leakage_mask_a = sf::st_union(roads_leakage, slope_leakage_mask)
leakage_mask_b = sf::st_union(leakage_mask_a, wetlands_leakage_mask)
leakage_mask_c = sf::st_union(leakage_mask_b, protected_areas_leakage_mask)

# save
sf::st_write(leakage_mask_c, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Leakage Masks/", delete_dsn=T)

Tally Features

road_count_whole = sf::st_intersection(road, leakage_belt_whole)
road_count = sf::st_intersection(road, leakage_belt)
road_length_whole = sum(sf::st_length(road_count_whole)) + sum(sf::st_length(road_count_whole))
road_length = sum(sf::st_length(road_count)) + sum(sf::st_length(road_count))
road_length_whole
road_length

waterways_count_whole = sf::st_intersection(waterways, leakage_belt_whole)
waterways_count = sf::st_intersection(waterways, leakage_belt)
waterways_length_whole = sum(sf::st_length(waterways_count_whole))
waterways_length = sum(sf::st_length(waterways_count))
waterways_length_whole
waterways_length

places_count_whole = sf::st_intersection(places, leakage_belt_whole)
places_count = sf::st_intersection(places, leakage_belt)
places_count_whole
places_count